#include <iostream>
#include <fstream>
#include <hash_map>
#include <vector>
#include <string>
#include <sstream>
#include <utility>
#include <ctime>
#include <iterator>
#include <algorithm>

using namespace std;

typedef int variable_type;
typedef hash_map<variable_type,double> joint_type;

const int ymask = 0xFFE00000;
const int xmask = 0x001FFFFF;
const int yshift = 32-10-1;

//const long long ymask = 0xFFE0000000000000;
//const long long xmask = 0x001FFFFFFFFFFFFF;


inline double mutual_information(const vector<vector<int> >& data, 
	const vector<int>& target, 
	const vector<int>& test_vector)

{
	hash_map<string, double> xmap, ymap, xymap;

	size_t N = target.size();
	size_t test_per_row = test_vector.size();
	ostringstream xstr, ystr, xystr;
	for(int i = 0; i < N; ++i)
	{
		for(int j = 0; j < test_per_row; ++j)
		{
			xstr << data[i][test_vector[j]];
		}
		ystr << target[i];
		xystr << xstr.str() << "," << ystr.str();

		pair<hash_map<string,double>::iterator, bool> ret = xmap.insert(make_pair(xstr.str(), double(1/double(N))));
		if(!ret.second)
		{
			ret.first->second += double(1/double(N));
		}

		ret = ymap.insert(make_pair(ystr.str(), double(1/double(N))));
		if(!ret.second)
		{
			ret.first->second += double(1/double(N));
		}

		ret = xymap.insert(make_pair(xystr.str(), double(1/double(N))));
		if(!ret.second)
		{
			ret.first->second += double(1/double(N));
		}

		xstr.str("");
		ystr.str("");
		xystr.str("");
	}

	double mi = 0.f;

	for(hash_map<string, double>::const_iterator it = xymap.begin();
		it != xymap.end(); ++it)
	{
		string xystring(it->first);
		size_t sep_pos = xystring.find_first_of(',');
		string xstring(xystring.substr(0, sep_pos));
		string ystring(xystring.substr(sep_pos+1));

		double px = xmap[xstring];
		double py = ymap[ystring];
		double pxy = it->second;

		mi += pxy*log(pxy/(px*py));
	}

	return mi;
}

inline double mutual_information_int(const vector<vector<int> >& data, 
	const vector<int>& target, 
	const vector<int>& test_vector)

{
	joint_type xmap, ymap, xymap;

	size_t N = target.size();
	size_t test_per_row = test_vector.size();
	variable_type xstr = 0, ystr = 0, xystr = 0;
	variable_type tmp = 0;
	for(int i = 0; i < N; ++i)
	{
		for(int j = 0; j < test_per_row; ++j)
		{
			tmp = data[i][test_vector[j]];
			tmp <<= (j*3);
			xstr |= tmp;
		}
		ystr = target[i];
		ystr <<= yshift;
		xystr = (ystr | xstr);

		pair<joint_type::iterator, bool> ret = xmap.insert(make_pair(xstr, 1.f));//double(1/double(N))
		if(!ret.second)
		{
			ret.first->second += 1.f;//double(1/double(N));
		}

		ret = ymap.insert(make_pair(ystr, 1.f));
		if(!ret.second)
		{
			ret.first->second += 1.f;//double(1/double(N));
		}

		ret = xymap.insert(make_pair(xystr, 1.f));
		if(!ret.second)
		{
			ret.first->second += 1.f;//double(1/double(N));
		}

		xstr = 0;
		ystr = 0;
		xystr = 0;
	}

	double mi = 0.f;

	for(joint_type::const_iterator it = xymap.begin();
		it != xymap.end(); ++it)
	{
		xystr = (it->first);
		ystr = (xystr&ymask);
		xstr = (xystr&xmask);

		double px = xmap[xstr]/double(N);
		double py = ymap[ystr]/double(N);
		double pxy = it->second/double(N);

		mi += pxy*log(pxy/(px*py));
	}

	return mi;
}


inline double mutual_information_sort(const vector<vector<int> >& data, 
	const vector<int>& target, 
	const vector<int>& test_vector)

{
	size_t N = target.size();
	size_t test_per_row = test_vector.size();
	variable_type xstr = 0, ystr = 0, xystr = 0;
	variable_type tmp = 0;
	vector<variable_type> x_in_data(N),y_in_data(N),xy_in_data(N);

	for(int i = 0; i < N; ++i)
	{
		for(int j = 0; j < test_per_row; ++j)
		{
			tmp = data[i][test_vector[j]];
			tmp <<= (j*3);
			xstr |= tmp;
		}
		ystr = target[i];
		ystr <<= yshift;
		xystr = (ystr | xstr);

		x_in_data[i] = xstr;
		y_in_data[i] = ystr;
		xy_in_data[i] = xystr;

		xstr = 0;
		ystr = 0;
		xystr = 0;
	}

	vector<double> x_prob_vec, y_prob_vec, xy_prob_vec;
	vector<variable_type> x_u, y_u, xy_u;
	x_u.reserve(N);
	y_u.reserve(N);
	xy_u.reserve(N);

	x_prob_vec.reserve(N);
	y_prob_vec.reserve(N);
	xy_prob_vec.reserve(N);

	sort(x_in_data.begin(), x_in_data.end());
	sort(y_in_data.begin(), y_in_data.end());
	sort(xy_in_data.begin(), xy_in_data.end());

	xstr = x_in_data[0];
	ystr = y_in_data[0];
	xystr = xy_in_data[0];

	x_u.push_back(xstr);
	y_u.push_back(ystr);
	xy_u.push_back(xystr);
	int xcount = 1, ycount = 1, xycount = 1;
	for(int i = 1; i < N; ++i)
	{
		if(xstr!=x_in_data[i])
		{
			xstr = x_in_data[i];
			x_u.push_back(xstr);
			x_prob_vec.push_back(double(xcount)/double(N));
			xcount = 0;
		}
		++xcount;

		if(ystr!=y_in_data[i])
		{
			ystr = y_in_data[i];
			y_u.push_back(ystr);
			y_prob_vec.push_back(double(ycount)/double(N));
			ycount = 0;
		}
		++ycount;

		if(xystr!=xy_in_data[i])
		{
			xystr = xy_in_data[i];
			xy_u.push_back(xystr);
			xy_prob_vec.push_back(double(xycount)/double(N));
			xycount = 0;
		}
		++xycount;
	}

	// last one
	x_prob_vec.push_back(double(xcount)/double(N));
	y_prob_vec.push_back(double(ycount)/double(N));
	xy_prob_vec.push_back(double(xycount)/double(N));

	double mi = 0.f;

	size_t num_xy = xy_u.size();
	vector<variable_type>::const_iterator x_u_begin = x_u.begin();
	vector<variable_type>::const_iterator x_u_end = x_u.end();
	vector<variable_type>::const_iterator y_u_begin = y_u.begin();
	vector<variable_type>::const_iterator y_u_end = y_u.end();
	for(size_t i = 0; i < num_xy; ++i)
	{
		vector<variable_type>::const_iterator it = lower_bound(x_u_begin, x_u_end, xy_u[i]&xmask);
		double px = x_prob_vec[distance(x_u_begin, it)];
		it = lower_bound(y_u_begin, y_u_end, xy_u[i]&ymask);
		double py = y_prob_vec[distance(y_u_begin, it)];

		double pxy = xy_prob_vec[i];
		mi += pxy*log(pxy/(px*py));
	}

	return mi;
}

vector<double> calculate_scores(const vector<vector<int> >& data, 
	const vector<int>& target, 
	const vector<vector<int> >& test_matrix)
{
	size_t num_tests = test_matrix.size();
	vector<double> scores(num_tests);

	for(size_t i = 0; i < num_tests; ++i)
	{
		scores[i] =  mutual_information_int(data, target, test_matrix[i]);
	}

	return scores;
}



vector<int> read_target_file(const string& file_name)
{
	ifstream input(file_name.c_str(), ios::binary);
	vector<int> target;
	if(!input)
	{
		cerr << "File " << file_name << " does not exist" << endl;
		return target;
	}

	target.reserve(10000);
	int num = 0;
	while(input>>num)
	{
		target.push_back(num);
	}

	return target;
}

vector<vector<int> > read_data_file(const string& file_name)
{
	ifstream input(file_name.c_str(), ios::binary);
	vector<vector<int> > data;
	if(!input)
	{
		cerr << "File " << file_name << " does not exist" << endl;
		return data;
	}

	data.reserve(1000);
	string line;
	while(getline(input, line))
	{
		int num;
		istringstream line_stream(line);
		vector<int> row;
		row.reserve(1000);
		while(line_stream>>num)
			row.push_back(num);

		data.push_back(row);
	}

	return data;
}


int main(int argc, char* argv[])
{
	srand(0);
	/*string target_file_name(argv[1]);
	string data_file_name(argv[2]);
	string matrix_file_name(argv[3]);
	string score_file_name(argv[4]);

	ofstream output_file(score_file_name.c_str(), ios::binary);
	*/
	int num_data = 1000, num_rows=2000;
	vector<int> target;// = read_target_file(target_file_name);
	vector<vector<int> > data;// = read_data_file(data_file_name);
	vector<vector<int> > matrix;// = read_data_file(matrix_file_name);
	for(int i = 0; i < num_data; ++i)
	{
		target.push_back(rand()%1024);
		data.push_back(vector<int>(500));
		for(int j = 0; j < 500; ++j)
		{
			data[i][j] = rand()%5;
		}
	}

	for(int i = 0; i < num_rows; ++i)
	{
		matrix.push_back(vector<int>(4));
		for(int j = 0; j < 4; ++j)
		{
			matrix[i][rand()%4] = rand()%500;
		}
	}

	//vector<int> target = read_target_file(target_file_name);
	//vector<vector<int> > data = read_data_file(data_file_name);
	//vector<vector<int> > matrix = read_data_file(matrix_file_name);
	clock_t t1,t2;
	t1=clock();
	vector<double> scores = calculate_scores(data, target, matrix);
	t2=clock();
	double diff ((double)t2-(double)t1);
	cout<<diff/CLOCKS_PER_SEC<<endl;



	ofstream output_file("new_method.txt", ios::binary);
	output_file.precision(5);	
	output_file.setf(ios::fixed);
	output_file.setf(ios::showpoint);
	copy(scores.begin(), scores.end(), ostream_iterator<double>(output_file, "\n"));
	return 0;
}